Introduction

Pour cette deuxième séance, nous allons apprendre à utiliser un planifateur d’itinéraire sur R afin de calculer des isochrones autour de chaque supermarché et ainsi pouvoir identifier les zones de chalandises.

1. Cartographier des isodistances

Une première information simple à obtenir mais très informative peut être obtenue en traçant des courbes d’isodistances autour de chaque supermarchés. On peut ainsi indentifier les zones non-deserviés.

Représentez les courbes d’isodistances (st_buffer(), tm_borders())

télécharger le package en .zip en cliquant sur ce lien
https://cran.r-project.org/src/contrib/Archive/opentripplanner/opentripplanner_0.4.0.tar.gz

Et installer Java 8 sur votre ordinateur:
https://www.java.com/fr/download/

library(dplyr)
library(ggplot2)
library(sf)
library(readr)
library(readxl)
library(tmap)
library(tmaptools)
library(here)

otp_available <- requireNamespace("opentripplanner", quietly = TRUE)
if (!otp_available) {
  otp_pkg <- here("Cours", "opentripplanner_0.4.0.tar.gz")
  if (file.exists(otp_pkg)) {
    try(install.packages(otp_pkg, repos = NULL, type = "source"), silent = TRUE)
    otp_available <- requireNamespace("opentripplanner", quietly = TRUE)
  }
}

## Warning in install.packages(otp_pkg, repos = NULL, type = "source"): 'lib =
## "/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library"' is not
## writable

if (otp_available) {
  library(opentripplanner)
} else {
  message("Package 'opentripplanner' absent: les sections OTP/isochrones seront ignorees.")
}

# Chargement des données locales (mêmes sources que test.Rmd)
iris_shp <- st_read(
  dsn = here("Cours"),
  layer = "georef-france-iris-millesime"
)

## Reading layer `georef-france-iris-millesime' from data source 
##   `/Users/loic/Documents/UGA/donnee_entreprise/ProjetECO/Cours' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 728 features and 30 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 4.74204 ymin: 44.69597 xmax: 6.359009 ymax: 45.88339
## Geodetic CRS:  WGS 84

sirene_shp <- st_read(
  dsn = here("Cours"),
  layer = "dataseed-sirene-2"
)

## Reading layer `dataseed-sirene-2' from data source 
##   `/Users/loic/Documents/UGA/donnee_entreprise/ProjetECO/Cours' 
##   using driver `ESRI Shapefile'
## replacing null geometries with empty geometries
## Simple feature collection with 342 features and 106 fields (with 5 geometries empty)
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 4.767436 ymin: 44.90643 xmax: 6.121839 ymax: 45.81497
## Geodetic CRS:  WGS 84

iris_data <- read_excel(here("Cours", "iris_isere.xlsx"))

# Zone d'étude: Grenoble + communes nord
sirene_coms <- c(
  "GRENOBLE", "MEYLAN", "LA TRONCHE",
  "CORENC", "SAINT-MARTIN-LE-VINOUX", "SAINT-EGREVE"
)
iris_coms <- c(
  "Grenoble", "Meylan", "La Tronche",
  "Corenc", "Saint-Martin-le-Vinoux", "Saint-Égrève"
)

sirene_grenoble <- sirene_shp %>%
  filter(libellecomm %in% sirene_coms & etatadminis == "Actif")
iris_grenoble <- iris_shp %>%
  filter(com_name %in% iris_coms)

# Préparation des variables IRIS (même logique que test.Rmd)
colnames(iris_grenoble)[20] <- "CODE_IRIS"

iris_data$nb_residence_princ <- iris_data$Res_princ_30_moins_40_m2 +
  iris_data$Res_princ_40_moins_60_m2 +
  iris_data$Res_princ_60_moins_80_m2 +
  iris_data$Res_princ_80_moins_100_m2 +
  iris_data$Res_princ_100_moins_120_m2

vars_to_keep <- c(
  "CODE_IRIS",
  "1er_quartile_euro", "3e_quartile_euro",
  "Actifs_15-64_ans2", "Pop_15-64_ans",
  "Delta_Pop_20-64_ans", "nb_residence_princ"
)

iris_grenoble <- merge(
  iris_grenoble,
  iris_data[, vars_to_keep],
  by = "CODE_IRIS",
  all.x = TRUE
)

iris_grenoble <- iris_grenoble %>%
  rename(
    premier_quartile_revenu = `1er_quartile_euro`,
    troisieme_quartile_revenu = `3e_quartile_euro`,
    Actif_total = `Actifs_15-64_ans2`,
    Pop_15_64 = `Pop_15-64_ans`,
    Delta_pop_20_64ans = `Delta_Pop_20-64_ans`,
    Nb_residence_principale = nb_residence_princ
  )

iris_grenoble <- iris_grenoble %>%
  mutate(
    across(
      c(
        premier_quartile_revenu,
        troisieme_quartile_revenu,
        Actif_total,
        Pop_15_64,
        Delta_pop_20_64ans,
        Nb_residence_principale
      ),
      ~ na_if(., 0)
    )
  )

st_crs(sirene_grenoble)

## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["latitude",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["longitude",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4326]]

# Buffer de 500 m autour des supermarchés
sirene_grenoble_l93 <- st_transform(sirene_grenoble, 2154)
buffer <- st_buffer(sirene_grenoble_l93, 500)
buffer <- st_transform(buffer, st_crs(iris_grenoble))

tmap_mode("view")
tm_shape(iris_grenoble) +
  tm_text(text = "iris_name", size = 0.8) +
  tm_polygons(
    "Pop_15_64",
    alpha = 0.5,
    style = "pretty",
    id = "iris_name",
    title = "Population (15-64 ans)"
  ) +
  tm_shape(buffer) +
  tm_borders("chartreuse3", lwd = 3) +
  tm_shape(sirene_grenoble) +
  tm_dots(
    "enseigne1et",
    legend.show = TRUE,
    title = "Enseigne",
    id = "enseigne1et",
    popup.vars = "trancheeffe.2",
    alpha = 0.3
  ) +
  tm_layout(legend.outside = TRUE)

Avec des cercles pleins plutôt que des courbes

tm_shape(iris_grenoble) +
  tm_text(text = "iris_name", size = 0.8) +
  tm_polygons(
    "Pop_15_64",
    alpha = 0.5,
    style = "pretty",
    id = "iris_name",
    title = "Population (15-64 ans)"
  ) +
  tm_shape(buffer) +
  tm_polygons("enseigne1et", legend.show = FALSE, id = "enseigne1et", alpha = 0.5) +
  tm_shape(sirene_grenoble) +
  tm_dots(
    "enseigne1et",
    legend.show = TRUE,
    title = "Enseigne",
    id = "enseigne1et",
    popup.vars = "trancheeffe.2",
    alpha = 0.8
  ) +
  tm_layout(legend.outside = TRUE)

##

## ── tmap v3 code detected ───────────────────────────────────────────────────────

## [v3->v4] `tm_polygons()`: instead of `style = "pretty"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style' to 'tm_scale_intervals(<HERE>)'
## [v3->v4] `tm_polygons()`: use `fill_alpha` instead of `alpha`.
## [v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_polygons()`: use `fill_alpha` instead of `alpha`.
## [v3->v4] `tm_polygons()`: use `fill.legend = tm_legend_hide()` instead of
## `legend.show = FALSE`.
## [v3->v4] `tm_dots()`: use `fill_alpha` instead of `alpha`.
## [tm_dots()] Arguments `legend.show` and `title` unknown.

2. Obtenir des isochrones

2.1 Mise en place du calculateur d’itinéraire

Nous allons utiliser OpenTripPlanner, un logiciel codé en Java et utilisé notamment dans le calculateur d’itinéraire du TAG. Pour fonctionner, OTP a besoin de 3 choses:

Un dossier & un sous-dossier specifiques à créer dans votre R working directory > OTP > graphs > default

Des données sur le réseau de transport GTFS (General Transit Feed Specification) à déposer dans OTP > graphs > default

Un fond de carte openstreetmap à déposer dans OTP > graphs > default

données GTFS

gtfs_img_candidates <- c(
  here("Cours", "gtfs.png"),
  here("Cours", "Capture.PNG"),
  here("Cours", "Capture.png")
)
gtfs_img <- gtfs_img_candidates[file.exists(gtfs_img_candidates)][1]
if (!is.na(gtfs_img)) {
  knitr::include_graphics(gtfs_img)
} else {
  cat("Ajoutez une capture GTFS dans Cours/gtfs.png pour l'affichage.")
}

Vous pouvez par exemple trouver les données GTFS sur le site M mobilité du TAG: data.mobilites-m.fr
Une fois téléchargé le fichier SEM-GTFS.zip renommez le en gtfs.zip et placez le dans OTP > graphs > default

fond de carte OSM

Il y a aussi beaucoup de sources, par exemple :
https://download.geofabrik.de/europe/france/rhone-alpes.html
Placez le fichier dans OTP > graphs > default

Mise en place d’OTP

path_data <- NULL
for (candidate in c(here("Cours", "OTP"), here("Cours", "data"))) {
  if (
    dir.exists(file.path(candidate, "graphs", "default")) ||
    file.exists(file.path(candidate, "Graph.obj"))
  ) {
    path_data <- candidate
    break
  }
}
if (is.null(path_data)) {
  path_data <- here("Cours", "OTP")
  dir.create(path_data, recursive = TRUE, showWarnings = FALSE)
}
path_otp <- otp_dl_jar()

# Vous n'avez besoin de lancer la construction du graphe qu'une fois
# log <- otp_build_graph(otp = path_otp, dir = path_data, memory = 15000)

log1 <- otp_setup(otp = path_otp, dir = path_data)
otpcon <- otp_connect()

otp_img_candidates <- c(
  here("Cours", "OTP_capture.png"),
  here("Cours", "OTP_capture.PNG")
)
otp_img <- otp_img_candidates[file.exists(otp_img_candidates)][1]
if (!is.na(otp_img)) {
  knitr::include_graphics(otp_img)
} else {
  cat("Ajoutez une capture OTP dans Cours/OTP_capture.png pour l'affichage.")
}

2.2 Estimer les isochrones (otp_isochrone())

iso <- otp_isochrone(
  otpcon,
  fromPlace = sirene_grenoble,
  fromID = sirene_grenoble$siret,
  mode = c("WALK", "TRANSIT"),
  maxWalkDistance = 2000,
  date_time = as.POSIXct(strptime("2022-02-22 08:35", "%Y-%m-%d %H:%M")),
  cutoffSec = c(5, 7, 9, 11) * 60
)
iso$minutes <- iso$time / 60

# Vérification
colnames(iso)[3] <- "siret"
iso <- merge(
  iso,
  st_drop_geometry(sirene_grenoble[, c("siret", "enseigne1et")]),
  by = "siret",
  all.y = TRUE
) %>%
  unique()

# Option de correction ponctuelle (si un point ne retourne pas d'isochrone)
if (nrow(sirene_grenoble) >= 16) {
  st_geometry(sirene_grenoble)[16] <- st_sfc(
    st_point(c(5.728253270386639, 45.17659752397455)),
    crs = st_crs(sirene_grenoble)
  )
}

# On recommence
iso <- otp_isochrone(
  otpcon,
  fromPlace = sirene_grenoble,
  fromID = sirene_grenoble$siret,
  mode = c("WALK", "TRANSIT"),
  maxWalkDistance = 1500,
  date_time = as.POSIXct(strptime("2022-02-22 08:35", "%Y-%m-%d %H:%M")),
  cutoffSec = c(5, 7, 9, 11) * 60
)
colnames(iso)[3] <- "siret"

iso <- merge(
  iso,
  st_drop_geometry(sirene_grenoble[, c("siret", "enseigne1et")]),
  by = "siret",
  all.y = TRUE
) %>%
  unique()
iso$minutes <- iso$time / 60

3. Cartes finales et indices

3.1 Carte intéractive avec isochrones

tmap_mode("view")
tm_shape(iris_grenoble) +
  tm_polygons(
    "Pop_15_64",
    alpha = 0.5,
    style = "pretty",
    id = "iris_name",
    title = "Population (15-64 ans)"
  ) +
  tm_shape(sirene_grenoble) +
  tm_dots(
    "enseigne1et",
    legend.show = TRUE,
    title = "Enseigne",
    size = 0.1,
    id = "enseigne1et",
    popup.vars = "trancheeffe.2"
  ) +
  tm_layout(legend.outside = TRUE) +
  tm_shape(iso) +
  tm_fill(
    "minutes",
    breaks = c(0, 5.01, 7.01, 9.01, 11.01),
    title = "Isochrones (minutes)",
    style = "fixed",
    labels = c("0 to 5", "5 to 7", "7 to 9", "9 to 11"),
    palette = "-BuPu",
    id = "minutes",
    alpha = 0.3
  ) +
  tm_borders()

3.2 Analyse de la concurrence

Calculez le nombre de supermarchés accessibles par IRIS (st_intersects())

iso420 <- iso %>% filter(time == 420)

iris_grenoble$nb_supermarche_per_iris <- rowSums(
  ifelse(st_intersects(x = iris_grenoble, y = iso420, sparse = FALSE), 1, 0)
)

tm_shape(iris_grenoble) +
  tm_polygons(
    "nb_supermarche_per_iris",
    style = "pretty",
    id = "iris_name",
    labels = c("0", "1", "2", "3", "4", "5"),
    title = "Nb de Supermarchés accessible en 7min"
  ) +
  tm_text(text = "iris_name", size = 0.8)

Bonus pour la 3ème session, on calcule le nombre de supermarchés de chaque enseigne accessibles par IRIS

data_temp <- matrix(0, nrow = nrow(iris_grenoble), ncol = nrow(iso420))
for (j in 1:nrow(iso420)) {
  for (i in 1:nrow(iris_grenoble)) {
    data_temp[i, j] <- as.character(st_intersects(iris_grenoble, iso420, sparse = FALSE)[i, j])
    data_temp[i, j] <- ifelse(data_temp[i, j] == "TRUE", iso420$enseigne1et[j], 0)
  }
}

data_temp2 <- matrix(0, nrow = nrow(data_temp), ncol = length(unique(iso420$enseigne1et)))
for (j in 1:length(unique(iso420$enseigne1et))) {
  data_temp2[, j] <- apply(
    data_temp,
    1,
    function(x) length(which(x == unique(iso420$enseigne1et)[j]))
  )
}

data_temp2 <- as.data.frame(data_temp2)
colnames(data_temp2) <- unique(iso420$enseigne1et)
write.csv(data_temp2, here("Cours", "nb_pdv_iris.csv"), row.names = FALSE)

Calculez la surface cumulée couverte par un supermarché pour chaque iris (st_intersection(), st_area())

intersect_pct <- st_intersection(iris_grenoble, iso420) %>%
  mutate(intersect_area = st_area(.)) %>%
  dplyr::select(iris_name, intersect_area) %>%
  group_by(iris_name) %>%
  mutate(intersect_area = sum(intersect_area)) %>%
  st_drop_geometry() %>%
  unique()

Puis calculez la surface de chaque IRIS, le taux de couverture par iris et faîtes une carte

iris_grenoble <- mutate(iris_grenoble, iris_area = st_area(iris_grenoble))

iris_grenoble <- merge(iris_grenoble, intersect_pct, by = "iris_name", all.x = TRUE)

iris_grenoble <- iris_grenoble %>%
  mutate(intersect_area = ifelse(is.na(intersect_area), 0, intersect_area)) %>%
  mutate(couverture = as.numeric(intersect_area / iris_area) * 100)

tm_shape(iris_grenoble) +
  tm_polygons(
    col = "couverture",
    alpha = 0.5,
    style = "cont",
    id = "iris_name",
    title = "Couverture (%)"
  ) +
  tm_text(text = "iris_name", size = 0.8) +
  tm_shape(sirene_grenoble) +
  tm_dots(
    "enseigne1et",
    legend.show = TRUE,
    title = "Enseigne",
    size = 0.1,
    id = "enseigne1et",
    popup.vars = "trancheeffe.2"
  ) +
  tm_layout(legend.outside = TRUE)

Enregistrer iris_grenoble pour la prochaine séance

st_write(
  iris_grenoble,
  here("Cours", "iris_grenoble2.gpkg"),
  delete_dsn = TRUE
)

Devoirs pour la prochaine séance

Préparez plusieurs cartes avec des isochrones différentes et d’autres variables au niveau des IRIS (carte 3.1) afin de proposer une localisation potentielle d’un nouveau supermarché à Grenoble. Proposez aussi une autre carte (3.2) construite à partir d’une isochrone différente (>7min). Ce travail sera évalué.

Geomarketing avec R - Session 2 (Grenoble + Zone Nord)